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Abstract 

A  two  dimensional  mesoscale  atmospheric  model  is  presented  based  on  two  systems 
of  dynamical  equations,  the  non-hydrostatic  compressible  Euler  equations  and  the 
corresponding  hydrostatic  system.  Both  equation  systems  are  discretized  with  a 
discontinuous  Galerkin  method  in  space  and  a  linear  semi-implicit  multistep  method 
in  time.  To  cover  elementwise  polynomial  spaces  up  to  order  four,  exact  quadrature 
rules  are  used.  The  method  is  applied  to  a  terrain  following  quadrilateral  grid.  To 
validate  the  models,  gravity  wave  propagation  and  mountain  wave  experiments  are 
analyzed.  The  numerical  models  exhibit  the  proper  wave  propagation  characteristics 
and  high  order  convergence  rates. 


1  Introduction 

The  compressible  Euler  equations,  or  the  non-hydrostatic  equation  system,  are  the 
governing  equations  for  atmospheric  motion  of  adiabatic  dry  and  frictionless  air. 
They  cover  important  dynamical  features  of  atmospheric  fluid  motion,  their  dynam¬ 
ics  is  extremely  non-linear  and  characterized  by  interactions  between  processes  of 
multiple  spatial  and  temporal  scales.  The  presence  of  sound  waves,  the  potential 
for  the  development  of  energetic  shocks,  and  the  potential  for  wave  breaking  and 
scale  collapses  are  fundamental  features  of  these  equations.  Sound  waves  are  the 
fastest  waves  contained  in  these  equations  and  usually  constitute  a  severe  CFL- 
condition  for  the  explicit  time  step  size  according  to  linear  stability  analysis.  Thus, 
filtering  the  dynamical  equations  to  exclude  sound  wave  propagation  is  one  way  to 
circumvent  this  problem.  E.g.  incompressible,  anelastic  and  pseudo-incompressible 
equations  do  not  contain  fast  sound  waves,  see  e.g.  Ogura  and  Phillips  (1962),  Dur- 
ran  (1989).  The  hydrostatic  approximation  leads  to  a  vertically  filtered  system,  but 
allows  horizontally  propagating  fast  Lamb  waves,  see  e.g.  Kalnay  (2003).  This  is 
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one  reason  why  for  several  decades,  the  hydrostatic  system  was  used  for  weather 
forecast  and  climate  simulations.  Another  way  to  avoid  strict  CFL-conditions  is  to 
consider  the  unfiltered  equation  system  and  to  use  implicit  or  semi-implicit  time 
integration  schemes  with  larger  stability  regions  appropriate  for  stiff  problems. 
Within  the  present  article,  we  consider  examples  for  both  types  of  equations.  These 
are  the  unfiltered  fully  non-hydrostatic  system  and  the  vertically  filtered  hydro¬ 
static  system.  For  these  two  systems  we  analyze  wave  propagation  characteristics 
caused  by  initial  perturbations  and  mountain  wave  forcings.  Specifically,  this  anal¬ 
ysis  estimates  the  range  of  horizontal  scales  that  is  properly  represented  within  the 
hydrostatic  equation  system.  It  should  be  noted  that  constructing  a  hydrostatic 
model  based  on  a  discontinuous  Galerkin  method  is  a  mathematically  challenging 
problem  that  we  show  here  how  to  solve  in  a  mathematically  consistent  approach. 
Further,  our  interest  is  to  present  a  detailed  convergence  study  to  compare  experi¬ 
mental  and  theoretical  spatial  convergence  rates.  To  our  knowledge,  this  has  never 
been  shown  before  for  either  hydrostatic  or  non- hydrostatic  atmospheric  models. 
For  2-dimensional  mesoscale  atmospheric  models,  a  number  of  experimental  setups 
are  known  from  the  literature.  Skamarock  and  Klemp  (1994)  studied  the  propaga¬ 
tion  of  inertia  gravity  waves  in  non- hydrostatic  and  hydrostatic  setups.  Meanwhile 
Pinty  et  al.  (1995)  give  examples  for  single  mountains  with  hydrostatic  and  non¬ 
hydrostatic  shape  size  and  Schar  et  al.  (2002)  propose  multiple  mountains  in  a 
consecutive  formation. 

In  our  current  paper,  we  discretize  the  equations  in  space  with  the  discontinuous 
Galerkin  method  (DGM),  one  possible  high  order  generalization  of  finite  volume 
methods.  Since  the  DGM  was  proposed  in  the  early  1970s  in  Reed  and  Hill  (1973), 
it  has  become  a  powerful  computational  tool  for  inviscid  and  viscous  problems  in 
gas  dynamics.  For  atmospheric  applications  Giraldo  et  al.  (2002),  Nair  et  al.  (2005), 
Giraldo  (2006),  Giraldo  and  Restelli  (2008),  Lauter  et  al.  (2008)  and  Restelli  and 
Giraldo  (2009)  have  successfully  applied  the  DGM.  It  features  high  order  accuracy, 
discrete  conservation  properties,  applicability  on  structured  and  unstructured  grids, 
and  robustness  especially  for  non-linear  problems.  The  locality  of  memory  usage 
offers  very  good  scalability  on  multi-core  architectures  and  on  graphics  processors, 
see  e.g.  Biswas  et  al.  (1994)  and  Kldckner  et  al.  (2009). 

Giraldo  and  Restelli  (2008)  and  Restelli  and  Giraldo  (2009)  have  presented  a  DGM 
solving  the  non-hydrostatic  system  in  a  mesocale  atmospheric  channel  on  a  rect¬ 
angular  grid.  The  grid  is  deformed  into  the  vertical  direction  to  obtain  a  terrain 
following  grid  structure.  The  temporal  evolution  of  the  conserved  variables  is  com¬ 
puted  solving  a  system  of  conservation  laws.  For  the  hydrostatic  equations,  the 
situation  is  different  because  the  vertical  momentum  component  is  no  longer  prog¬ 
nostic.  Further,  these  equations  consist  of  the  conservation  laws  for  the  remaining 
prognostic  variables  with  the  hydrostatic  constraint  to  determine  the  diagnostic  ver¬ 
tical  momentum  IF  as  a  Lagrangian  multiplier.  For  the  conservation  laws,  the  DGM 
is  applied  using  the  standard  Rusanov  numerical  flux.  For  the  hydrostatic  equations 
this  numerical  flux  is  modified  considering  jump  contributions  of  IF  only  into  the 
vertical  direction.  Our  approach  shows  that  the  non-hydrostatic  and  hydrostatic 
equations  can  be  implemented  in  a  unified  way  and  their  differences  are  controlled 
by  a  hydrostatic  switch  parameter  5h-  As  described  above  semi-implicit  (implicit- 
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explicit)  time  discretization  schemes  offer  improved  stability  properties  and  lead  to 
larger  time  steps  even  for  stiff  systems  containing  fast  waves.  For  the  DGM  Dolc- 
jsi  and  Feistauer  (2004),  Kanevsky  et  al.  (2007),  Giraldo  and  Restelli  (2009)  and 
Restelli  and  Giraldo  (2009)  have  shown  the  successful  application  of  semi-implicit 
Runge-Kutta  methods  and  linear  multi-step  methods.  In  our  current  work  for  both 
systems,  the  non-hydrostatic  and  hydrostatic,  a  2nd  order  semi-implicit  linear  multi- 
step  method  is  chosen  that  combines  robust  stability  and  high  order  properties  in 
time. 

The  article  is  structured  the  following  way.  After  the  introduction  of  the  non¬ 
hydrostatic  and  hydrostatic  systems  in  section  2,  section  3  describes  the  DGM  and 
the  temporal  scheme  for  both  equation  systems.  In  section  4  the  experimental  results 
including  a  convergence  study  are  presented.  Finally,  in  section  5  we  summarize  the 
main  results  and  propose  future  work. 

2  Non-Hydrostatic  and  Hydrostatic  Systems 

In  this  section  the  governing  dynamical  equations  are  discussed  representing  the 
mathematical  model  describing  the  atmospheric  flow  problems  considered  in  this 
article.  Two  different  systems  of  equations  for  compressible  dry  inviscid  air  with 
gravitational  acceleration  will  be  considered,  the  non-hydrostatic  and  hydrostatic 
equation  system. 

Because  our  approach  is  to  apply  a  discontinuous  Galerkin  method  to  the  equations, 
we  restrict  our  attention  to  equations  in  conservation/flux  form.  Two  flux  formu¬ 
lations  of  the  non-hydrostatic  system  are  known:  one  uses  density,  momentum  and 
potential  temperature  (the  ©-formulation)  and  the  other  uses  density,  momentum 
and  total  energy  as  prognostic  variables.  For  the  hydrostatic  approximation,  the 
prognostic  equation  for  the  vertical  momentum  W  is  modified  into  a  diagnostic  one. 
As  a  consequence,  the  kinetic  and  total  energy  equations  lose  their  prognostic  struc¬ 
ture  whereas  the  potential  temperature  equation  remains  prognostic.  Because  of 
this  observation,  we  restrict  ourselves  to  the  ©-formulation  of  the  non- hydrostatic 
system. 

The  governing  equations  will  be  considered  in  the  2-dimensional  spatial  domain 

D  =  {(a;,  z)  G  (. xL ,  xR)  x  R  |  zB{x)  <  z  <  zT} 

where  Xl  <  Xr  are  the  x-components  of  the  left  and  right  boundary,  zr  '■  ( Xl ,  Xr)  — > 
M  is  a  given  height  function  describing  the  domain  bottom  orography  and  Zt  is  the 
z-component  of  the  rigid  lid  on  top.  In  Q,  we  will  denote  coordinates  by  (aq,^) 
and  ( x ,  z)  depending  on  the  context. 

2.1  Physical  Variables 

The  non-hydrostatic  system  and  its  hydrostatic  approximation  are  given  by 

A Hdtq  +  div  f (q)  =  r(q )  in  D  x  M>0  (1) 
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with  the  conserved  variables  q  =  (p,U,W,Q)T  =  (p,  pu,  pw,  p9)T ,  density  p,  the 
velocity  vector  ( u,w ),  the  potential  temperature  9  and  the  time  domain  —  {t  e 
M  1 1  >  0}.  The  flux  function  and  the  right  hand  side  are  defined  by 
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with  the  pressure  p  =  c o©7,  Co  =  (. R/Pq )7,  the  gas  constant  i?,  a  constant  pressure 
Po,  the  isentropic  exponent  7,  the  Poisson  constant  k  and  the  gravitational  constant 
g.  z)  >  0  is  the  prescribed  Raleigh  damping  parameter  and  pa,  Ua,  @CT  the  cor¬ 
responding  fields  that  realize  the  non-reflecting  boundary  within  a  sponge  layer,  see 
section  4.  For  the  parameter  switch  5h  =  1,  Eq.  (1)  describes  the  non- hydrostatic 
system,  for  5h  =  0  these  are  the  equations  in  hydrostatic  approximation. 


2.2  Perturbation  Variables 


To  enhance  the  hydrostatic  equilibrium  of  the  model  the  numerical  method  will  be 
applied  to  the  perturbation  variable  with  respect  to  a  reference  field.  For  that  we 
assume  a  hydrostatic  reference  field  qr  =  (pr  =  pa,  0, 0,  ©r  =  ©o-)T,  pr  =  Co©7, 
9r  =  y-  with  the  properties  dzpr  =  — gpr ,  dtqr  =  0,  dxqr  =  0.  The  choice  for  the 
reference  fields  depends  on  the  numerical  experiment  and  will  be  reported  in  section 
4.  Defining  the  perturbation  variables 


q'  =  q~  Qr,  P\z,  q')  =  p(qr(z)  +  q')  -  pr{z ), 
system  (1)  can  be  formulated  in  terms  of  q'  by 

A  Hdtq'  +  div  f\z,  q)  =  r(q')  in  D  x  M>0 


with 
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To  simplify  the  notation  for  the  rest  of  the  article,  we  omit  the  primes. 


/  -V  \ 
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(2) 


3  Conservative  Numerical  Method 

To  approximate  the  spatial  derivatives  of  the  dynamical  equations  (2),  we  consider 
a  discontinuous  Galerkin  method  on  a  logically  rectangular  grid  combined  with  a 
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second  order  semi-implicit  backward  difference  time  step  method.  The  grid  geometry 
is  defined  by  a  terrain  following  ay-coordinate  giving  a  deformation  into  the  vertical 
direction.  Using  the  concept  of  isoparametric  finite  elements,  the  discrete  function 
space  consists  of  high  order  tensor  product  polynomials  in  each  element.  Based  on 
an  integral  form  of  the  equation  system,  a  high  order  discontinuous  Galerkin  method 
is  defined.  Thus,  the  conservation  properties  for  the  conserved  variables  are  assigned 
to  the  discrete  variables. 

3.1  Terrain  following  Grid 

The  numerical  method  is  based  on  a  rectangular  computational  grid  based  on  a 
terrain  following  ay-coordinate.  This  vertical  coordinate  with  a  linear  decay  of 
terrain  slopes,  see  e.g.  Gal-Chen  and  Somerville  (1975),  is  defined  by  a  global 
coordinate  mapping 

7  :  O  — >■  Q,  7(x,  z)  —  (  x,  zB(x)— — -  +  x 

\  Zt 

from  the  computational  domain  O  =  (xl,  xr)  x  (0,  zt)  to  the  spatial  domain  Q.  For 
the  horizonal  and  vertical  element  numbers  Nx  and  Nz,  a  rectangular  grid  structure 
is  defined  on  O  by  the  supporting  points 


xL  =  x0  <  Xi  <  ..  <  xNx  =  xR,  0  =  z0  <  z1  <  ..  <  zNz  =  zT. 

For  each  index  (i,j)  =  {1, ..,  Nx}  x  {1, ..,  Nz}  the  computational  and  its  correspond¬ 
ing  spatial  element  are  defined  by 

E c,ij  (*U— 1  j  Xi)  X  (77_i ,  7; ) ,  Eij  7 (-^c,p), 

such  that  O  is  covered  by  all  elements 

n=  ij  e 

EeT 

of  the  grids  triangulation  T  =  {Eij}i=it_^Nxj=i^,tNz- 

3.2  Discontinuous  Galerkin  Method 

On  each  element  E  G  T  with  respect  to  the  reference  element  Ec  C  O,  the  isopara¬ 
metric  tensor  product  polynomial  space  of  degrees  at  most  k  >  1  is  defined  by 

Q\E)  =  :  E  R  |  ip  o  7|Ec  e  Qkc(Ec)}, 

Qc(Ec )  =  span {Pij  :  Ec  R  \  0  <  i,  j  <  k]j)t]{x,  z)  =  xlz3}. 

We  notice  the  dimension  of  the  polynomial  space  to  be  Nk  =  dim(5A'(i7)  =  (k  +  l)2 
and  the  set  of  integer  numbers  Nk  =  {1  ,..,Nk}.  Thus,  the  coordinate  mapping  7 
defines  both  the  geometry  of  the  element  E  and  the  polynomial  space  on  E.  Based 
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on  the  polynomial  space  Qk(E),  the  discrete  discontinuous  function  space  is  defined 
by 

V  =  {veL°°(n)\v\EeQk(E),  WE  G  T}. 

Because  p  G  V  are  polynomials  on  each  grid  element  E  G  T,  the  condition  99  G  L°° 
does  not  constitute  an  additional  constraint  to  the  discrete  functions. 

The  starting  point  for  the  space-discrete  formulation  is  an  appropriate  integral  form 
of  the  dynamic  equations  (2).  This  is  obtained  multiplying  (2)  with  a  smooth 
(continuous  in  hi  with  derivatives)  test  function  <f>  :  — >  M4,  assuming  a  smooth 

solution  q  of  (2),  integrating  over  each  E  G  T  and  applying  integration  by  parts, 
that  is 

E  dE  E 

Assuming  the  indices  i  —  1,2  and  j  =  1,  ..,4,  the  terms  are  defined  by  /  :  V<F  = 
Yj1)3  fijdxj^i,  $  '  (/Fe)  =  $ifijVE,j  and  vE  G  M2  is  the  normal  vector  outward 
on  dE.  The  discrete  version  of  this  integral  form  gives  the  discontinuous  Galerkin 
method,  the  condition  that  the  space-discrete  solution  q(.,t)  G  V4  has  to  fulfill, 
namely 

(AHdtq,$)L2m  +  F(q(.,t),$)  =  0  for  all  $Gh4  (3) 

with  the  discontinuous  Galerkin  operator  F(q,&)  =  J^EeT  Fe(q,  for  divf(q)  — 
r(q )  and 


FE(q,  <I>)  -  ->.7  /  :  V$  +  $  •  r  dx  +  J  •  f(z ,  qin ,  qout,  vE)  da.  (4) 

E  dE 

Because  the  values  of  q  are  not  uniquely  defined  on  dE,  inner  and  outer  values 
qin  and  qout  are  considered  and  the  flux  function  is  substituted  by  a  numerical  flux 
function  /.  For  the  boundary  integral  along  dE,  at  the  element  boundary  point 
xq  G  dE,  the  inner  value  is  defined  by 

qin(x  0)  =  Ihn  q(x). 

x£E,x^x  0 

For  an  inner  element  edge  point  Xq  G  dE  \  d£l  respectively  a  boundary  edge  point 
x0  G  dE  D  <9f2,  the  outer  value  is  defined  by 


qout(x  0)  =  lim  q(x)  resp.  qout(x  0) 

x€.£l\E  ,x— >xq 


( 


lim 

x€E,x^x  0 
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\ 

•  V  V 


This  way,  reflecting  boundary  conditions  on  dtt  are  prescribed  for  the  flow.  For  the 
boundary  integral  in  (4),  a  numerical  flux  function  replaces  the  flux  f(q).  Among 
the  different  choices  for  the  numerical  flux,  we  have  taken  the  Rusanov  numerical 
flux 

-  1 

f(z ,  qi,  <?2,  u)  =  2  l/(z’  qi">u  +  ~  +  Id4  ^)(g2  -  qi)]  (5) 
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where  we  have  used  the  decomposition  v  =  (i/x,  vz)  and  Id4  is  the  rank-4  identity 
matrix.  Because  for  the  hydrostatic  case  5h  =  0  the  vertical  momentum  jumps 
are  weighted  with  z/?,  this  term  vanishes  along  vertically  aligned  element  faces. 
To  represent  the  space-discrete  solution  q(.,t)  G  V4  of  the  discontinuous  Galerkin 
method  (3),  a  basis  of  V  is  specified.  For  that,  in  each  element  E  E  T,  a  Lagrangian 
basis  (tpE,i)ieNk  is  defined  with  respect  to  a  tensor  product  Gauss- Lobatto  points 
as  collocation  points,  see  Bos  et  al.  (2000).  Let  us  continue  each  basis  function 
(fE,i  :  E  — >  M  by  zero  to  the  whole  domain  and  denote  this  discontinuous  function 
again  by  (pE,i-  Then,  (< PE,i)EeT,ieNk  establishes  a  Lagrangian  basis  of  V.  Now,  each 
scalar  component  qi  of  the  discrete  solution  q  can  be  represented  with  respect  to 
this  basis,  that  is 

Nk 

qi(x,  z,t)  =  qi,E,i(t)(pE,i(x,  z ), 

1=1 

with  the  components  qitE,i  and  the  component  vector  q  =  {qi,E,i)EeT-i=i ..  4-ieNk- 
Applying  this  to  (3)  yields  the  quasilinear  differential  algebraic  system 


M”  ft  +  =  °- 


(M 

0 

0 

\0 


0  0 
M  0 
0  8hM 
0  0 


°\ 

0 

0 

M) 


(6) 


with  the  mass  matrix  M  =  (fn<PE,i<PF,j) (E,i),(F,j)eTxNk  and  the  appropriate  vector 
function  F .  For  the  non- hydrostatic  system  (5h  =  1),  (6)  simplifies  to  a  system  of 
ordinary  differential  equations  by  inverting  the  mass  matrix  Mh- 
The  construction  of  the  vector  F  includes  FE(q,  <£>)  in  (4)  and  further  the  evaluation 
of  integrals  over  E  and  dE.  This  is  done  applying  appropriate  quadrature  rules. 
For  that,  we  follow  Cockburn  and  Shu  (2001)  to  obtain  k  +  1-order  formal  accuracy 
and  consider  tensor  Gauss-Lobatto  quadrature  rules  of  order  2k  in  each  element  E. 
On  each  edge  of  E,  Gauss-Lobatto  rules  of  order  2k  +  1  are  applied  that  imply  the 
evaluation  of  the  integrand  at  the  quadrature  points,  and  differ  from  the  collocations 
points  of  the  Lagrange  polynomial  basis. 


3.3  Semi-Implicit  Linear  Multistep  Method 

One  possibility  to  discretize  the  semi-discrete  problem  (6)  is  a  linear  multistep 
method.  Due  to  advection  terms  and  the  non-linear  wave  speed  A  in  the  Rusanov 
flux,  the  discontinuous  Galerkin  operator  F  is  non-linear.  Although  Newton-type 
solvers  for  this  non-linear  problem  are  a  possible  choice,  numerical  efficiency  is  still 
an  open  problem.  Here,  a  semi-implicit  (or  implicit-explicit)  multistep  method  is 
chosen,  to  take  advantage  of  a  large  stability  region,  see  Hundsdorfer  and  Verwer 
(2003),  Giraldo  (2005),  Giraldo  et  al.  (2009).  Further,  the  hydrostatic  constraint 
can  be  enforced  implicitly  in  every  time  step;  this  (along  with  the  flux  function 
defined  in  (5))  is  the  key  to  constructing  a  mathematically-consistent  discontinuous 
Galerkin  method  for  the  hydrostatic  equations. 

To  apply  the  semi-implicit  method,  a  linear  discontinuous  Galerkin  operator  L  is  in¬ 
troduced  to  approximate  F.  In  section  3.2,  the  non-linear  operator  F  is  constructed 
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based  on  the  flux  function  /.  Now,  the  same  procedure  is  chosen  to  construct  L 
based  on  a  linearized  flux  function  /.  With  the  hydrostatic  reference  fields  pr,  9r , 
©r  in  section  2.2,  a  Taylor  approximation  of  /  yields 


l(z,q) 


/  u 

w  \ 

Pi 

0 

0 

pi 

\eru 

erw) 

Pi(z,  ©)  =  c07©/  1(^)0. 


With  respect  to  l  a  linear  discontinuous  Galerkin  operator  for  div/(g)  —  r(q)  is 
defined  by 


£(?,*)  =  £- 
EeT 


l  :  V$  +  $  •  r  dx  +  /  <I>m  •  l(z,  qin,  qout ,  uE)  da 


dE 


with  a  constant  speed  of  sound  Ao  >  0  and  the  linear  Rusanov  numerical  flux 
-  1 

l(z,  qu  q2,  v)  =  2  [Z(a  +  l(z’  ~  Ao(A Hv2x  +  Id4  ^)(?2  -  qi)]  ■ 

Different  to  the  Rusanov  flux  in  (5),  the  wave  speed  A0  is  a  constant  and  ensures 
the  linearity  of  the  numerical  flux  l.  Finally,  the  linear  operator  L  is  constructed 
from  L(q,  <f>)  applying  the  representation  of  the  Lagrangian  basis  of  V.  Now,  the 
differential  algebraic  system  (6)  can  be  rewritten  with  a  linear  L  and  a  non-linear 
F  —  L,  that  is 

+  L(q)  +  (F-  L)(q)  =  0. 

For  the  temporal  discretization,  the  second  order  semi-implicit  backward  difference 
formula  BDF2,  a  linear  3-step  method,  is  used  treating  L  implicitly  and  F  —  L 
explicitly,  see  Giraldo  and  Restclli  (2009).  For  each  time  step  tn  — >  tn+1  this  can  be 
written  as 


Mh  (  r+1  -  OLm<r~m )  +  qAt  J2  Pm(F  -  L)(qn~m )  +  VA tL(<T+1)  =  0  (7) 

\  m= 0  /  m= 0 

with  the  parameters  a0  =  |,  a.\  —  —  rj  —  |,  j30  —  2  and  (3\  =  —1.  Giraldo 
and  Restelli  (2009)  have  further  studied  a  set  of  higher  order  implicit  BDFs,  which 
would  be  desirable  for  an  high  order  application.  But,  because  the  stability  regions 
of  the  higher  order  BDFs  do  exclude  a  neighborhood  of  0  on  the  imaginary  axis, 
these  alternatives  tend  to  generate  instabilities. 

The  linear  operator  L  is  chosen,  such  that  the  implicit  operator  L  includes  the 
fast /sound  waves  and  explicit  operator  F  —  L  the  advective  waves.  The  A-stability  of 
the  implicit  part  of  the  BDF  circumvents  a  time  step  restriction.  But  to  obtain  linear 
stability,  the  explicit  part  of  the  BDF  method  causes  a  maximum  CFL-timestep 
A tad  with  respect  to  the  advective  wave  speed.  Because  in  section  4  we  concentrate 
on  accuracy  rather  than  on  stability,  the  time  steps  used  for  the  experiments  are 
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far  below  the  stability  limit  Atad-  To  control  the  model  time  step,  we  define  a 
characteristic  time  that  takes  the  sound  waves  to  pass  a  grid  element,  that  is 


A  tso  = 


Ax 


cso(2fc  +  1) 

with  the  sound  wave  speed  cso  =  y/u2  +  w2  +  \fypfp.  The  motivation  for  this 
choice  is  the  CFL-condition  in  Cockburn  and  Shu  (2001)  for  explicit  Runge-Kutta 
Discontinuous-Galerkin  methods  both  of  order  k  +  1  applied  to  the  one-dimensional 
linear  case. 


3.4  Discrete  Hydrostatic  Balance 


Both  options  Sh  —  0  and  6h  =  1  have  been  considered  with  a  backward  differ¬ 
encing  formula  in  the  last  section.  For  the  non-hydrostatic  case,  Eq.  (7)  recovers 
the  approach  of  Restelli  and  Giraldo  (2009)  but  for  the  energy  form  of  the  Euler 
equations. 

The  hydrostatic  system  case  with  8h  —  0  has  a  different  structure.  While,  the 
equations  for  p,  U  and  0  in  (7)  remain  the  same,  the  third  equation  in  (7)  does  not 
include  a  discrete  time  derivative  for  d \W .  How  does  this  equation  look  like?  To 
answer  this  question,  we  consider  the  test  function  <f>  =  (0,  0,  <p,  0)  and  observe 


(F  -  L)(q,  $)  =  ^  -  (p  -  pi)dz(f  dx  +  /  ipin(p  -  Pi)vE,zda, 
EeT  i  L 


£(9,*)  =  E-  / 

E&T  I 


dE 

Pidz(f-  gpp  dx  +  J  <PinPiVE,z  da 

E  BE 

where  p  and  pi  are  the  values  of  the  pressure  functions  on  the  element  boundary 
defined  by  the  Rusanov  numerical  fluxes.  Applying  this  to  (7)  for  all  E  G  T  and 
<p  G  V  yields 


-  Pm(P  ~  +  Pl+1)9^  dx 

g  m= 0 

+  /  (^2  /3m(p  ~  Pi)n~m  +  P?+1)vvE,zda  =  -  j  gpn+1ip  dx. 
BE  m=0  E 

Provided  the  temporal  truncation  error  can  be  neglected 

l 

p"+1  =  5>„(p-p,)”_m+Pr+1  +  0(Ai2), 

m= 0 

the  third  equation  in  (7)  can  be  considered  as  a  discrete  hydrostatic  balance  at  time 
tn+1,  namely 


J  pn+1dz(p  dx  + 


J  pn+1(pisE,z  da 


J  gpn+1<fdx. 


(9) 


E  BE  E 

In  other  words,  Eq.  (9)  is  the  conservative  discontinuous  Galerkin  representation  of 
the  hydrostatic  balance. 
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4  Model  Validation 


The  experimental  results  presented  within  this  section  will  show  the  applicability  of 
the  discontinuous  Galerkin  method  to  the  non-hydrostatic  and  hydrostatic  systems. 
Characteristic  wave  features  can  be  established  and  convergence  studies  show  the 
expected  convergence  rates  for  smooth  numerical  solutions.  The  experiments  are 
performed  for  the  polynomial  orders  k  =  1,2,  3, 4,  with  the  linear  semi-implicit 
multistep  method  BDF2. 

For  each  model  experiment,  the  model  resolution  is  given  by  a  number  hp  >  0. 
Depending  on  the  anisotropy  parameter  [3  >  0,  hp  reflects  different  resolutions  into 
the  horizontal  and  vertical  direction.  To  specify  this,  the  anisotropic  element  size  of 
E  €  T,  the  anisotropic  grid  size  and  the  anisotropic  model  resolution  are  defined  by 

\E\p  —  max(  max  \x$  —  x\\,  (3  max  \zo~zi\), 

( xo,zo),{xi,zi)gE  (xo,zo),(xi,zi)eE 

a  |  .  Axp 

Axp  =  max\E\p,  hp  = 

with  the  polynomial  order  k.  This  way,  a  grid  with  h  =  10km  and  [3  =  2  effectively 
has  a  horizontal  resolution  of  10km  and  a  vertical  resolution  of  5km.  We  will  write 
h  instead  of  hp  if  the  assignment  is  clear. 


4.1  Gravity  Waves 


The  first  validation  test  is  a  gravity  wave  caused  by  an  initial  potential  temperature 
perturbation  O'  of  an  uniformly  stratified  atmosphere.  Following  Giraldo  and  Restelli 
(2008)  for  the  potential  temperature  9,  the  initial  condition  is  9r  +  6'  with  the 
reference  field  9r(z)  =  90exp(N2  z  /  g) ,  the  Brunt-Vaisala  frequency  N  =  0.01  s^1, 
the  temperature  constant  9q  =  300K  and  the  horizontal  velocity  u  =  20m/s.  The 
initial  perturbation  is  given  by 

9\x,  z)  =  A 60  sin(/^)4(a,  xc,  x)  (10) 


with  A 90  =  0.01K,  l 


7r/10km,  the  Agnesi  ’’versiera”  mountain  profile 


A(a,  xc,  x)  = 


a2  +  (x  —  xc )2  ’ 


the  center  xc  and  the  width  parameter  a.  In  a  linearized  Boussinesq  atmosphere, 
Skamarock  and  Klemp  (1994)  have  derived  an  analytic  solution  for  the  potential 
temperature  perturbation 


9un(x,  z,t)  =  sm(lz)A90a  /  exp(— ka)  cos(Af)  cos(xk)  dk 


(ii) 


with  A  =  Nk/^Snk2  +  P  that  gives  a  qualitative  approximation  to  the  solutions  of 
the  non-linear  non- hydrostatic  system  as  well  as  for  the  hydrostatic  system  in  (1). 
Because  of  the  dispersion  relation 


c(k) 


N 

VSHk2  + 1 2 


(12) 
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for  the  Boussinesq  atmosphere,  we  expect  the  fastest  waves  in  the  non-hydrostatic 
system  case  5h  —  1  to  have  the  same  speed  N /l  as  the  dispersion  free  waves  in  the 
hydrostatic  system  Sh  =  0. 

This  test  is  used  for  two  setups,  in  short  and  long  atmospheric  channels.  In  the 
short  channel  setup,  we  consider  the  spatial  domain  hi  =  (0,300km)  x  (0,10km), 
the  constants  xc  =  90km,  a  =  5km  and  the  integration  time  T  =  3000s.  In  the 
long  channel  setup,  the  domain  =  (0,  6000km)  x  (0, 10km)  with  the  constants 
xc  =  1800km,  a  =  100km  and  the  integration  time  T  =  60000s  is  used.  For  all 
experiments  with  the  short  channel  setup,  the  anisotropic  model  resolution  is  defined 
by  the  anisotropy  parameter  (3  =  1.  For  the  long  channel  setup  the  parameter  is  set 
to  /3  =  20. 

For  the  non- hydrostatic  and  hydrostatic  Boussinesq  equations,  Fig.  1  depicts  the  po¬ 
tential  temperature  perturbation  after  3000s  for  the  short  channel  setup  and  after 
60000s  for  the  long  channel  setup.  For  the  long  channel,  because  the  differences 
between  the  non-hydrostatic  and  hydrostatic  cases  are  not  visible,  the  hydrostatic 
plot  is  omitted.  The  numerical  solutions  of  (2),  for  the  potential  temperature  per¬ 
turbation  9  with  k  —  1,  h  —  0.6km,  At  =  A tso  are  given  in  Fig.  2.  This  setup  yields 
a  time  step  of  At  =  1.1s.  Again  for  the  long  channel,  no  differences  are  visible 
between  the  non-hydrostatic  and  hydrostatic  systems  and  we  only  plot  the  first  one. 
Qualitatively,  the  experiments  for  the  non-hydrostatic  and  hydrostatic  systems  are 
in  good  agreement  with  the  Boussinesq  solutions  in  Fig.  1;  nonetheless,  visible  dif¬ 
ferences  can  be  observed,  especially  in  the  vertical  structure.  Further,  the  linear 
dispersion  property  (12)  is  in  very  good  agreement  with  the  non-linear  experiments. 
This  can  be  seen  by  the  dispersion  for  the  short  channel  and  the  non- hydrostatic 
system  with  8h  —  1.  For  the  long  channel  setup,  the  perturbation  is  dominated  by 
long  horizontal  wave  numbers  and  thus  the  wave  dispersion  effect  is  reduced.  For  the 
hydrostatic  system  in  both  channels  the  experiments  show  no  significant  dispersion, 
as  expected. 

For  smooth  solutions,  the  convergence  properties  of  the  discontinuous  Galerkin 
method  are  studied.  Given  a  fixed  polynomial  order  k,  the  order  of  convergence 
k  +  1  for  the  L2-norm  can  be  expected.  To  generate  smooth  solutions,  the  initial 
condition  (10)  has  to  be  modified.  First,  the  periodic  boundary  conditions  lead  to 
jumps  along  the  lateral  boundaries  of  fh  Second,  the  reflecting  boundary  conditions 
on  top  and  bottom  of  can  be  reinterpreted  as  a  vertical  continuation  by  mirror¬ 
ing  (10)  along  the  horizontal  axis  that  features  jumps  in  the  vertical  derivatives  of 
(10).  Thus  only  for  the  convergence  study,  (10)  is  replaced  by  the  smooth  potential 
temperature  perturbation 


0'sm(x,  z)  =  A6>0(1  -  cos(2 lz))A{a,  xc,  x)<pXc(x  -  xc) 

with  (pe(x)  =  exp( — tz^z)  if  M  <  |s|  and  <pe(x)  =  0  otherwise. 

The  relative  L2-error  rj  for  the  numerical  solution  9  is  derived  with  respect  to  a 
reference  solution  6Q  with  the  model  resolution  h,  that  is 


V  =  V{h,  0) 


\\0  - 

|| Or  +  @h\  \  L2(fl) 
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Figure  1:  Gravity  waves,  section  4.1,  potential  temperature  perturbation  in  (11)  for 
Boussinesq  equations,  contour  interval  5  x  10-4K,  top:  short  channel,  T  =  3000s, 
non-hydrostatic  system,  middle:  short  channel,  T  =  3000s,  hydrostatic  system, 
bottom:  long  channel,  T  =  60000s,  non-hydrostatic  system;  In  the  long  channel,  the 
non-hydrostatic  and  hydrostatic  plots  look  the  same.  Thus,  the  hydrostatic  plot  is 
omitted. 
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x  [km] 


Figure  2:  Gravity  waves,  section  4.1,  model  results  for  potential  temperature  pertur¬ 
bation  6 ,  contour  interval  5  x  10-4K,  top:  short  channel,  T  =  3000s,  non-hydrostatic 
system,  middle:  short  channel,  T  =  3000s,  hydrostatic  system,  bottom:  long  chan¬ 
nel,  T  =  60000s,  non- hydrostatic  system;  In  the  long  channel,  the  non- hydrostatic 
and  hydrostatic  plots  look  the  same.  Thus,  the  hydrostatic  plot  is  omitted. 
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Figure  3:  Gravity  waves,  section  4.1,  L2-errors  and  convergence  rates  for  the  non¬ 
hydrostatic  system,  T  =  300s,  BDF2,  At  =  |A tso,  left:  short  channel,  right:  long 
channel. 


With  the  reference  value  9r  from  section  2.2  and  the  perturbation  value  9g,  9r  +  9h  is 
the  physical  value  for  the  potential  temperture.  Then,  the  convergence  rate  ag(h,  k) 
for  the  variable  9  is  defined  depending  on  the  model  resolution  h  and  the  fixed 
polynomial  order  k.  For  that,  the  numerical  solutions  9h. ,  92h,  9ih  are  considered 
with  respect  to  the  three  model  resolutions  h,  2 h  and  4 h.  Then 


a  =  otoih ,  k ) 


In  ri(h,  92h) 

ln(2h)  —  ln(4/i) 


gives  an  approximation  to  the  convergence  order  of  the  method. 

The  convergence  properties  of  the  spatial  discretization  are  analyzed  for  the  fixed 
integration  time  of  T  =  300s.  To  separate  the  temporal  error,  we  seek  the  largest 
time  step  size  for  the  BDF2  method  such  that  error  and  convergence  rates  do  not 
change  for  smaller  time  steps.  For  the  presented  experiments,  this  is  nearly  ob¬ 
tained  for  the  time  step  At  =  |A tso,  with  A tso  in  (8).  The  relative  L2-errors  r/ 
and  convergence  rates  a  are  illustrated  in  Fig.  3  for  the  non- hydrostatic  system  in 
both  channels.  The  model  converges,  that  is  the  error  decreases  monotonously  with 
increasing  model  resolution  and  increasing  polynomial  order.  Furthermore  for  each 
polynomial  order  k  —  1,  ..,4,  a  is  close  to  the  expected  value  k  +  1.  Only  for  the 
high  order  case  k  =  4,  a  is  slightly  limited  by  the  temporal  error.  In  contrast  to 
this,  Fig.  4  shows  the  same  errors  rj  and  rates  a  for  non-smooth  initial  conditions 
(10)  and  for  longer  time  steps,  respectively.  The  initial  conditions  cause  a  slight 
degradation  but  the  long  time  step  strongly  degrades  the  accuracy. 

The  same  convergence  experiments  but  for  the  hydrostatic  system  are  depicted  in 
Fig.  5.  Again,  the  method  converges  and  the  convergence  rates  match  very  well  with 
the  expected  values.  Compared  to  the  non- hydrostatic  system,  the  convergence  rates 
have  even  improved. 


4.2  Hydrostatic  Mountain  Waves 

Mountain  waves  arise  in  a  steady-state  flow  perturbed  by  a  mountain  orography. 
This  test  case  is  comprised  of  a  horizontal  wind  field  u  with  a  single  hydrostatic 
mountain  with  a  large  horizontal  size  a  such  that  the  forcing  frequency  u/a  is  small 
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Figure  4:  Gravity  waves,  section  4.1,  L2-errors  and  convergence  rates  for  the  non¬ 
hydrostatic  system,  short  channel,  T  =  300s,  BDF2,  left:  non-smooth  initial  condi¬ 
tion  (10),  At  =  | A tso,  right:  smooth  initial  condition,  At  =  5 A tso. 


Figure  5:  Gravity  waves,  section  4.1,  L2-errors  and  convergence  rates  for  the  hy¬ 
drostatic  system,  T  =  300s,  BDF2,  At  =  |A tso,  left:  short  channel,  right:  long 
channel. 
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Figure  6:  Hydrostatic  mountain,  section  4.2,  horizontal  velocity  perturbation  u, 
contour  interval  0.005m/s,  cutout  of  fl,  T  =  lOh,  left:  non- hydrostatic  system, 
right:  hydrostatic  system. 


compared  to  the  Brunt- Vaisala  frequency  N  =  g^/n/RT.  Thus,  internal  waves  are 
excited  and  propagate  freely  into  the  vertical  direction,  see  Smith  (1979).  We  use 
the  same  setup  as  in  Giraldo  and  Restelli  (2008)  which  is  similar  to  that  in  Pinty 
et  al.  (1995)  and  Klcmp  and  Lilly  (1978).  In  the  spatial  domain  Q  =  (0,240km)  x 
(0,  20km),  an  Agnesi  mountain  is  given  by 

h(x,  z)  =  hoA(a,  xc,  x ) 

with  ho  =  lm,  xc  =  120km  and  the  mountain  width  parameter  a  =  10km.  The 
initial  condition  is  a  hydrostatic  isothermal  atmosphere  with  T  =  250K,  p(z  =  0)  = 
lOOkPa  and  the  horizontal  velocity  u  =  20m/s.  The  hydrostatic  characteristics  of 
the  mountain  is  obtained  by  the  relation  aN/2nu  1.  On  the  bottom,  reflecting 
boundary  conditions  are  assumed  while  on  the  top  and  lateral  boundaries  sponge 
layers  of  6km  respectively  20km  thickness  are  applied.  The  lateral  sponges  enforce 
the  inflow  and  outflow  conditions  to  agree  with  the  initial  data. 

For  both  equation  systems,  the  solution  converges  to  a  quasi-stationary  state.  After 
a  simulation  time  of  T  —  lOh,  Fig.  6  depicts  the  horizontal  velocity  perturbations 
for  k  —  1,  h  —  1.2km,  (3  =  4.8  and  At  =  7.6s  in  a  cutout  of  Q.  For  the  non¬ 
hydrostatic  system,  the  wave  structure  agrees  very  well  with  Restelli  and  Giraldo 
(2009).  Further,  the  results  for  the  large  horizontal  wave  numbers  of  the  mountain 
profile  lead  to  a  very  good  agreement  between  the  non-hydrostatic  and  hydrostatic 
systems. 

The  convergence  analysis  for  rj  and  a  is  given  in  Fig.  7  for  the  integration  time 
T  =  360s  and  the  anisotropy  parameter  /l  =  1.5.  Convergence  is  observed  for 
the  non-hydrostatic  and  hydrostatic  systems,  but  only  for  resolutions  better  than 
approximately  6km.  The  convergence  rates  are  close  to  the  expected  values  but  the 
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Figure  7:  Hydrostatic  mountain,  section  4.2,  L2-errors  and  convergence  rates,  T  = 
360s,  BDF2,  At  =  A tso,  left:  non-hydrostatic  system,  right:  hydrostatic  system. 


higher  order  convergence  rates  (k  =  3, 4)  are  reduced  slightly.  This  effect  might  be 
attributed  to  the  non-smooth  acoustic  wave  fronts  during  the  transient  phase  forced 
by  the  mountain. 

4.3  Non- Hydrostatic  Mountain  Waves 

Non- hydrostatic  mountain  waves  are  excited  over  a  mountain  orography  with  a 
smaller  horizontal  size  a  than  in  section  4.2.  Here,  the  forcing  frequency  u/a  is  large 
compared  to  the  Brunt-Vaisala  frequency  N  such  that  the  vertical  wave  propagation 
is  damped  by  buoyancy  forces.  We  study  a  similar  test  setup  as  in  Giraldo  and 
Restclli  (2008).  In  the  spatial  domain  Q  =  (0, 100km)  x  (0,20km),  the  same  Agnesi 
mountain  profile  as  in  section  4.2  is  used  but  with  the  parameters  ho  =  lm,  xc  = 
50km  and  a  =  1km.  The  initial  condition  is  an  uniformly  stratified  atmosphere  with 
N  =  0.01s_1,  T(z  =  0)  =  280K,  p(z  =  0)  =  lOOkPa  and  the  horizontal  velocity 
u  =  lOm/s.  Now,  the  dimensionless  parameter  aN j2i\u  <C  1  describes  the  non¬ 
hydrostatic  characteristics  of  the  mountain.  The  same  boundary  conditions  as  in 
section  4.2  are  assumed,  but  with  sponge  layer  thicknesses  of  6km  (top)  and  10km 
(lateral).  For  both  equation  systems,  after  a  simulation  time  of  T  —  5h,  a  quasi¬ 
stationary  state  is  achieved.  Fig.  8  depicts  the  horizontal  velocity  perturbations  for 
k  —  1,  h  —  0.4km,  (3  =  2  and  At  =  4.2s  in  a  cutout  of  Q.  The  wave  structure  for  the 
non-hydrostatic  system  is  in  very  good  agreement  to  Restelli  and  Giraldo  (2009). 
For  the  hydrostatic  system,  the  amplified  horizontal  gradients  lead  to  strong  waves 
located  above  the  orographic  perturbation.  Thus,  the  small  horizontal  scales  of  the 
orographic  perturbation  yield  significant  differences  between  the  non- hydrostatic 
and  hydrostatic  systems. 

The  convergence  analysis  for  r)  and  a  is  given  in  Fig.  9  for  the  integration  time 
T  =  180s  and  the  anisotropy  parameter  (3  =  1.  Because  the  horizontal  mountain 
scale  is  10  times  smaller  compared  to  section  4.2,  convergence  is  established  for 
the  non-hydrostatic  and  hydrostatic  systems  below  a  model  resolutions  of  1km. 
The  convergence  rates  are  satisfactory,  but  below  the  expected  theoretical  rates  of 
convergence.  One  reason  is  the  non-smoothness  of  the  solution,  as  described  in 
section  4.2.  Another  reason  could  be  the  limited  spatial  resolution  which  reduces 
the  confidence  in  the  values  of  the  estimated  convergence  rates  a. 
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Figure  8:  Non-hydrostatic  mountain,  section  4.3,  horizontal  velocity  perturbation  u, 
contour  interval  O.OOlm/s,  cutout  of  il,  T  —  5h,  left:  non-hydrostatic  system,  right: 
hydrostatic  system. 


Figure  9:  Non-hydrostatic  mountain,  section  4.3,  L2-errors  and  convergence  rates, 
T  =  180s,  BDF2,  At  =  A tso,  left:  non- hydrostatic  system,  right:  hydrostatic  sys¬ 
tem. 
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This  section  and  section  4.2  allow  the  conclusion  that  the  hydrostatic  equation 
system  is  able  to  model  mountain  waves  very  good  above  the  horizontal  mountain 
scale  of  approximately  10km.  In  contrast  to  that,  finer  mountain  scales  lead  to 
artificial  wave  propagation.  One  possibility  to  prevent  this  error  in  a  hydrostatic 
model  is  to  apply  the  model  with  a  resolution  above  the  10km  scale.  Finer  model 
resolutions  are  appropriate  only  if  the  horizontal  mountain  scale  is  still  above  the 
10km  scale. 


4.4  Schar  Mountain  Waves 

An  experiment  constituting  a  superposition  of  hydrostatic  and  non-hydrostatic  moun¬ 
tain  waves  is  described  in  Schar  et  al.  (2002).  The  orography  is  given  by 


h(x)  =  h0ex  p 


(x  -  xcy 


cos 


7l(x  -  Xc) 

A 


with  the  mountain  height  ho  =  250m,  xc  =  25km  and  the  wave  parameters  a  =  5km 
and  A  =  4km.  In  the  spatial  domain  Q  —  (0,  50km)  x  (0,  20km)  for  the  initial  and 
reference  fields  an  uniformly  stratified  atmosphere  is  assumed  with  N  =  O.Ols^1, 
T(z  =  0)  =  280K,  p(z  =  0)  =  lOOkPa  and  the  horizontal  velocity  u  =  lOm/s.  The 
same  boundary  condition  as  in  the  previous  sections  are  used.  The  top  and  lateral 
sponge  layer  thicknesses  are  6km  and  5km,  respectively. 

For  the  non-hydrostatic  system  after  a  simulation  time  of  T  =  lOh,  the  quasi¬ 
stationary  horizontal  velocity  perturbations  are  plotted  in  Fig.  10  for  k  —  2,  h  = 
0.2km,  f3  =  0.5  and  At  =  1.9s.  Following  the  discussion  in  section  4.3  about  small 
scale  mountains  in  the  hydrostatic  system,  here  the  hydrostatic  results  are  omitted. 
Like  in  section  4.2  qualitatively  this  result  is  in  good  agreement  to  Schar  et  al. 
(2002)  and  Restclli  and  Giraldo  (2009).  Small  and  large  scale  waves  are  apparent, 
the  larger-scale  waves  propagate  and  the  smaller-scale  waves  rapidly  decay  with 
height. 

The  results  of  the  convergence  studies  for  the  non- hydrostatic  system  is  given  in 
the  right  column  of  Fig.  10  for  the  integration  time  T  =  360s  and  the  anisotropy 
parameter  [3  —  1.  Qualitatively  here,  the  convergence  results  are  similar  to  that 
in  section  4.3.  Convergence  is  obtained  below  a  model  resolution  of  1km.  More  so 
than  in  sections  4.2  and  4.3,  there  is  a  significant  difference  between  the  observed 
convergence  rates  and  their  expected  values. 


5  Summary 

The  article  presents  a  2-dimensional  mesoscale  model  of  the  atmosphere  based  on 
the  non- hydrostatic  system  and  its  hydrostatic  approximation.  Both  systems  are 
discretized  with  a  discontinuous  Galerkin  method  in  space  respecting  a  logically  rect¬ 
angular  terrain  following  grid.  We  present  an  unified  formulation  of  non-hydrostatic 
and  hydrostatic  atmospheric  systems  using  a  mathematically  consistent  discontinu¬ 
ous  Galerkin  formulation;  the  introduction  of  a  specialized  flux  function  is  one  of  two 
keys  to  the  success  of  this  approach.  To  avoid  small  time  steps  respecting  a  severe 
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Figure  10:  Schar  Mountain,  Section  4.4,  non-hydrostatic  system,  left:  horizontal 
velocity  perturbation  u,  contour  interval  0.2m/s,  T  =  lOh,  right:  L2-errors  and 
convergence  rates,  T  =  360s,  BDF2,  At  =  A tso. 

CFL-condition  a  semi-implicit  linear  multi-step  method  in  time  has  been  applied. 
The  time  step  size  is  controlled  according  to  a  heuristic  CFL-condition  respecting 
the  grid  size  and  the  polynomial  order  of  the  method.  The  second  key  concept 
in  constructing  an  unified  non- hydrostatic /hydrostatic  discontinuous  Galerkin  relies 
on  a  semi-implicit  time-integration  method  whereby  the  hydrostatic  constraint  is 
imposed  implicitly  as  part  of  the  time-integration  strategy. 

The  validation  experiments  show  correct  wave  propagation  properties  as  well  as 
the  expected  convergence  rates.  The  gravity  wave  propagation  experiments  show 
good  agreement  with  the  analytic  solution  of  the  linearized  Boussinesq  atmosphere. 
Specifically,  wave  dispersion  for  the  non-hydrostatic  system  and  no  dispersion  for 
the  hydrostatic  system  are  observed.  With  smooth  initial  conditions,  high  order 
convergence  rates  Axk+1  are  observed  for  polynomial  orders  k  =  1, ...,  4.  For  the 
non-hydrostatic  and  hydrostatic  systems,  the  results  for  the  linear  hydrostatic  moun¬ 
tain  waves  qualitatively  agree.  For  the  non- hydrostatic  system,  the  result  for  the 
linear  non- hydrostatic  mountain  waves  are  qualitatively  correct.  As  expected,  qual¬ 
itatively  the  results  of  the  hydrostatic  system  are  different.  Because  the  horizontal 
size  of  the  Schar  mountains  are  in  the  non-hydrostatic  scale,  differences  between  both 
systems  are  expected  and  are  indeed  shown  to  differ.  The  non- hydrostatic  system 
gives  the  well  known  wave  propagation  structure.  Although  the  experiments  show 
convergence  for  all  mountain  wave  experiments,  the  convergence  rates  are  below  the 
expected  theoretical  order.  Only,  the  non-hydrostatic  system  for  the  hydrostatic 
mountain  test  gives  the  convergence  order  Axfc+1  for  k  —  1,2. 

For  2-dimensional  physical  problems,  like  the  study  of  mountain  wave  propagation, 
the  presented  model  is  already  an  appropriate  tool.  But  for  more  realistic  mesoscale 
applications,  a  3-dimensional  model  is  required.  This  could  be  achieved  generalizing 
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the  present  code  to  3-dimensional  grids  based  on  cubic  or  prismatic  grid  elements. 
This  work  has  already  begun  and  we  shall  report  our  findings  for  more  compli¬ 
cated  flow  problems  in  the  future,  perhaps  including  moist  physics,  and  other  basisc 
physical  parameterizations. 
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